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results  obtained  by  this  (direct)  method  are  compared  to  estimates 
derived  from  applying  the  exponential  magnitude -frequency  relation- 
ship of  natural  seismicity.  It  is  found  that  the  latter  (indirect)  method 
in  general  gives  significantly  lower  90  percent  detection  thresholds. 
This  is  explained  theoretically  as  resulting  from  the  fact  that  the  in- 
direct method  fails  to  take  into  account  the  global  variance  of  observed 
signal  amplitudes  for  given  seismic  events.  It  is  concluded  that,  if  a 
good  reference  system  is  available,  the  direct  method  will  produce 
estimates  that  are  more  reliable  than  those  provided  by  the  indirect 
method. 
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ABSTRACT 


This  report  deals  with  the  general  problem  of  estimating  the 
incremental  detection  capability  of  a seismic  station,  i.  e.  , the  probability  of 
detection  as  a function  of  event  magnitude.  A maximum  likelihood  procedure 
is  introduced  to  estimate  the  station  detection  thresholds  by  comparison  with 
an  independent  reference  station  or  network.  Approximate  confidence  limits 
for  the  estimated  parameters  are  computed.  The  results  obtained  by  this 
(direct)  method  are  compared  to  estimates  derived  from  applying  the  expon- 
ential magnitude-frequency  relationship  of  natural  seismicity.  It  is  found 
that  the  latter  (indirect)  method  in  general  gives  significantly  lower  90  per- 
cent detection  thresholds.  This  is  explained  theoretically  as  resulting  from 
the  fact  that  the  indirect  meth  id  fails  to  take  into  account  the  global  variance 
of  observed  signal  amplitudes  for  given  seismic  events.  It  is  concluded  that, 
if  a good  reference  system  is  available,  the  direct  method  will  produce  esti- 
mates that  are  more  reliable  than  those  provided  by  the  indirect  method. 
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Technical  Applications  Center  will  be  responsible  for  inf<  rmation  contained 
herein  which  has  been  supplied  by  other  organizations  or  contractors,  and  this 
document  is  subject  to  later  revision  as  may  be  necessary.  The  views  and  con- 
clusions presented  are  those  of  the  authors  and  should  not  be  interpreted  as  ^ 
necessarily  representing  the  official  policies,  either  expressed  or  implied,  of 
the  Advanced  Research  Projects  Agency,  the  Air  Force  Technical  Applications 
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SECTION  1 


INTRODUCTION 


The  detection  capability  of  a seismic  station  or  network  for 
events  from  a specific  region  is  usually  referred  to  in  terms  of  its  incremen- 
tal detection  probability.  This  is  defined  as  the  probability  of  detecting  an 
event,  given  the  event  magnitude.  In  particular,  the  90  percent  detection 
threshold  is  often  quoted  as  a measure  of  performance;  this  is  the  magnitude 
at  which  the  station  is  expected  to  detect  90  percent  of  ail  events. 

Several  methods  have  been  devised  to  estimate  the  detection 
probability  function  of  a seismic  system.  In  general,  such  methods  can  be 
assigned  to  one  of  three  main  classes: 

• Estimates  based  on  seismic  noise  studies  - By  measuring 
the  seismic  noise  level,  estimating  the  signal-to-noise  ratio 
required  for  detection  and  assuming  a signal  variance,  one 
can  reasonably  well  predict  the  actual  detection  performance 
of  a system. 

• Estimates  based  on  seismicity  and  observed  detection  perfor- 
mance - This  is  a two  step  procedure.  First  the  seismicity 
of  a region  is  estimated  by  extrapolating  the  observed  data, 
using  the  exponential  magnitude-frequency  relationship.  Then 
the  observed  number  of  events  is  compared  to  the  estimated 
seismicity  in  order  to  establish  detection  thresholds. 

• Estimates  based  on  comparison  to  a reference  system  - A set 
of  events  reported  by  an  independent  reference  system  is  first 
selected.  The  percentage  actually  detected  at  each  magnitude 


r 
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by  the  station  in  question  is  then  used  to  obtain  threshold  esti 
mates. 


We  will  i.  this  report  refer  to  the  second  and  third  estimation 
mtthods  described  above  as  indirect  and  direct  estimation,  respectively. 

The  main  topic  of  this  report  is  to  present  a new  approach  to 
the  direct  method  of  estimation,  using  a maximum  likelihood  technique.  Ex- 
amples of  application  are  included,  and  the  results  are  compared  to  those  ob 
tained  by  other  methods. 

Section  II  of  this  report  establishes  a model  of  the  detection 
probability  function  which  has  been  found  useful  for  threshold  estimation.  In 
this  model,  the  probability  of  detecting  an  event  of  a given  magnitude  m is 
assumed  to  be  a cumulative  Gaussian  distribution  function: 


P(Detect  m) 


(1-D 


where  /t  and  <r  are  unknown  parameters.  It  is  shown  that  the  parameters 
should  be  interpreted  differently  according  to  which  method  of  estimation  ic 
being  used.  This  has  the  very  important  implication  that  different  methods 
of  estimation  may  be  expected  to  produce  different  results,  so  that  a careful 
in  erpretation  is  necessary. 

In  Section  III  the  likelihood  function  for  the  direct  estimation 
method  is  developed,  and  ah  iroximate  confidence  limits  for  the  estimated 
parameters  are  computed.  The  validity  of  the  approximations  is  examined 
by  applying  a simulation  model.  We  also  include  a brief  description  of  a 
maximum  likelihood  method  for  the  indirect  estimation  problem,  as  develop- 
ed by  Lacoss  and  Kelly  (1969).  We  choose  an  approach  which  is  slightly  dif- 
ferent from  theirs  in  order  to  show  that  no  hypothesis  of  Poisson  distribution 
of  natural  seismicity  is  required  to  develop  the  likelihood  function. 
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Section  IV  presents  some  examples  showing  how  the  direct  es- 
timation technique  can  be  applied  in  practical  situations.  A comparison  be- 
tween the  results  obtained  by  the  direct  and  indirect  methods  is  carried  out 
for  two  aftershock  sequences  recorded  by  the  NORSAR  and  LASA  arrays.  It 
is  found  that  the  two  methods  yield  significantly  different  results,  but  that  the 
apparent  disagreement  can  be  adequately  explained  by  the  considerations  pre- 
sented in  Section  II. 

main  results  from  this  report  are  summarized  in  Section 
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SECTION  II 

THE  GAUSSIAN  MODEL  FOR  EVENT  DETECTION  PROBABILITY 
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This  section  addresses  the  general  problem  of  establishing  a 
model  for  the  detection  probability  curve  of  a seismic  s/stem,  i.  e.  , the  pro- 
bability that  the  system  detects  an  event  as  a function  of  event  magnitude  m. 

The  case  of  single  station  or  seismic  array  detection  is  dis- 
cussed first,  and  it  is  shown  that  under  reasonable  assumptions  the  detection 
curve  will  have  the  form  of  a Gaussian  distribution  function. 

x 

2 

2-1/2  f 2 (T  . ...  , . 

P(  Detect/  m=x)  = (Zna  ) • I e dt  (II-l) 

-00 

This  holds  true  whether  the  magnitude  m is  defined  as  the  magni'ude  esti- 
mated by  a specific  station  or  as  a "true"  event  magnitude  measured  by  a 
hypothetical  "perfect"  network.  However,  the  resulting  values  of  the  para- 
meters M and  rr  will  be  different  in  the  two  cases. 

The  detection  probability  curve  of  a seismic  network  is  also 
discussed  briefly.  It  is  concluded  that  the  Gaussian  model  does  not  apply 
directly  to  this  case,  but  that  it  may  still  be  useful  as  an  approximation  to 
part  of  the  detection  curve. 

A.  SINGLE  STATION 

1.  Derivation  of  the  Detection  Probability 

In  order  to  determine  the  probability  of  a seismic  station  or 
array  detecting  an  event  from  a specified  region,  we  make  the  following  as- 
sumptions (illustrated  in  Figure  II-l): 


II-l 
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Parameters 


FIGURE  II- 1 

ILLUSTRATION  OF  THE  GAUSSIAN  HYPOTHESIS  OF  STATION  DETECTION 
THRESHOLD  MAGNITUDE  AND  STATION  EVENT  MAGNITUDE  DISTRIBUTION 


Given  that  an  event  has  true  magnitude  m=x  , then  the  magni- 


tude m of  the  corresponding  signal  arriving  at  the  station  is 

**  2 

a normal  variable:  ~ N(x  , o-^) 


The  event  Is  detected  at  the  station  provided  mR>mrj.  , where 


m is  a threshold  magnitude  determined  by  the  seismic  noise 


level  and  the  characteristics  of  t'  - detection  algorithm.  We 

2 

assume  that  m is  a normal  variable:  m ~ N(fi  ,w  ) . 


The  random  variables  m and  m are  assumed  independent. 

R T 


Under  these  assumptions,  the  probability  of  detection,  given 


m=x,  is  as  follows: 


P(  Detect/  m=x)  = P{m  >mT/m=x) 


= P(m  - m >0/m=x) 
K I 


(U-2) 


Since  the  conditional  distribution  of  (m^-m^),  given  m-x,  is  Gaussian,  we 


obtain: 


P(  Detect/ m=x)  = cfc 


(H-3) 


where 


2 2 2 
<r  = * <rR 


(H-4) 


and  <t  denotes  the  standard  cumulative  Gaussian  distribution  function. 


The  validity  of  this  model  for  station  detection  capability  rest  s, 
or  course,  with  the  validity  of  assuming  a normal  distribution  for  and 


For  the  station  magnitude  , the  normality  assertion  means 


that  for  a given  event,  world  wide  observed  magnitudes  follow  a 
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normal  distribution,  with  a standard  deviation  ^ R that  arises 
mainly  from  source  effects  and  differences  in  radiation  patterns. 
This  has  been  substantiated  by  several  studies  (Lanz,  1966; 
Freedman  1967).  Furthermore,  we  assume  that  the  variance 
of  this  distribution  is  independent  of  event  magnitude.  To  our 
knowledge,  this  has  not  been  explicitly  demonstrated  by  any 
author,  although  data  from  Bungum  and  Husebye  (1974)  com- 
paring NORSAR  and  PDE  magnitudes  give  some  support  to  this 
assertion. 

• For  the  threshold  magnitude  m^  , the  main  influencing  random 

factor  is  the  time  variability  of  the  noise  level  for  the  station. 
Evidence  for  a lognormal  distribution  of  the  seismic  noise  am- 
plitudes for  a given  station  as  a function  of  time  has  been  found 
by  Gerlach  et  al.  , (1966)  and  Alsup  and  Becker  G973)  and  is 
supported  by  results  from  Barnard  and  Whitel.-.w  (1972).  Ob- 
viously, a lognormal  distribution  of  r.oise  amplitudes  implies 
a normal  distribution  of  noise  "magnitudes". 

Besides  thv.  above  considerations,  we  mention  briefly  some 
additional  factors  that  influence  the  distributions  of  rr  R and  rv^.  Meas  *re- 
ment  errors  of  amplitude  and  period  are  obvious  examples,  another  potential 
source  is  systematic  errors  in  the  distance  factor.  The  spectral  characteris- 
tics of  the  received  signal  are  important  for  detectability  purposes,  and  here- 
fore  influence  n>y«  whether  an  automatic  detector  or  visual  inspection  is  ap- 
plied. In  the  automatic  detector  case,  our  normality  assumption  means  that 
filter  SNR  gain  (in  dB)  is  normally  distributed.  Similarly,  the  distribution  of 
beamforming  loss  influences  the  variance  of  m^  in  the  case  of  an  array  sta- 
tion. 

One  final  important  consideration  is  the  effects  of  applying  the 
Gaussian  model  to  events  within  a large  region.  In  this  case,  th<  normality 
assumption  is  only  valid  if  the  geographic  seismicity  distribution  of  the  region 
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is  such  that  the  observed  B-factors  are  normally  distributed.  This  is  clearly 
not  valid  in  general;  therefore,  it  is  necessary  to  require  that  the  model  only 
be  applied  to  a region  for  which  the  variation  in  B-factors  are  negligible  com- 
pared to  the  variation  of  the  other  components  in  our  random  model. 

This  last  requirement  is  not  very  restrictive,  for  example, 
over  the  epicentral  distance  30-80  degrees,  the  bodywave  distance  factor 
v aries  only  by  + 0.  1 mfa  units.  This  is  well  below  the  value  of  = 0.  4 
found  by  Veith  and  Clawson  (19V’.)  as  a typical  value  of  the  standard  deviation 
of  magnitude  measurements  at  single  sensor  stations,  assuming  that  standard 
B-factor  tables  are  used.  In  the  case  of  large  aperture  arrays,  such  as  LASA 
and  NORSAR,  one  can  expect  a lower  value  of  it  , due  to  the  fact  that  local 
effects  are  smoothed  out  when  the  signal  is  averaged  over  a number  of  sen- 
sors. A value  of  tr  betwee.i  0.25  and  0.30  seems  reasonable  in  these  two 

R 

cases;  thus  the  B-factor  variation  is  still  relatively  insignificant. 

2.  Examples  of  Application 

In  the  following  subsection  several  exampies  are  given  of  pos- 
sible situations  to  which  the  Gaussian  model  could  be  applied.  It  is  important 
to  note  the  conditional  nature  of  the  probability  expression  II-3,  since  this  is 
the  key  to  a proper  understanding  of  the  various  viewpoints  that  are  presented. 

To  be  specific,  let  us  analyze  three  cases;  dealing  with  the 
NORSAR  P-wave  detection  probability  in  terms  of  "true"  mb,  NORSAR  rr>b  and 
LASA  mb»  respectively.  We  will  assume  that: 

• Given  that  an  event  has  "true"  magnitude  m=x,  then 

the  distribution  of  the  NORSAR  magnitude  m is 

N,*  + V 'rN> 

the  diptr  ibution  of  the  LASA  magniutde  m is 

2 

N(x  + bL,  <tl) 


The  NORSAR  "threshold  magnitude"  mT  is  N(/iT,<7T) 

The  random  variables  m^,  rn^,  an<^  m-p  are  independent. 

a.  NORSAR  Detection  Curve  in  Terms  of  "True"  Magnit  ide  m 

2 

Given  that  m=x,  then  m.T  is  N(x+b  , tr  ),  so  that 

N N N 


P(  Detect/ m=x)  - P(mN>mT/  m=x) 


■*( 


(x  - MT)  + b 

(r, 


") 


(11-5) 


where 


2 , 2 

cr  + c* 

T N 


(11-6) 


b.  NORSAR  Detection  Curve  in  Terms  of  NORSAR  Magnitude 


mN‘ 


In  this  case,  the  problem  is  to  determine  the  probability  of  de- 
tection, given  that  a signal  arrives  at  NORSAR  corresponding  to  a magnitude 
= x.  (Thus  it  is  implicitly  assumed  that  the  signal  has  an  value  even 

if  it  is  not  detected.  ) 


Given  that  m = x , it  is  clear  that  detection  occurs  provided 
N 


x>m^,  , so  that 


/X'M 

P(  Detect/  m =x)  = <t>  I I 

N \ / 


(11-7) 


It  is  evident  that  an  indirect  procedure  is  required  in  order  to 
estimate  NORSAR  detectability  in  terms  of  m . One  such  procedure  would 

IN 

be  to  analyze  the  distribution  of  the  seismic  noise  level,  and  thereby  estimate 
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and  <r^.  Another  procedure  is  to  assume  that  the  number  of  earthquakes 
N(x)  exceeding  a given  NORSAR  magnitude  x may  be  expressed  by 


log  N(x)  = a - bx  (natural  logarithm). 


( II  - 8 ) 


and  use  this  expression  in  conjunction  with  the  observed  number  of  events  to 
determine  detection  thresholds,  (Lacoss  and  Kelly,  19f>9;  Bungum  and  Huse- 
bye,  1974). 

Note  that  the  formula  ( II- 8 ) is  more  commonly  expressed  in 
terms  of  base  10  logarithms.  If  this  is  done,  a and  b become  a'  = a/loglO 
and  b'  = b/loglO,  respectively. 

c.  NORSAR  Detection  Curve  in  Terms  of  LASA  Magnitude 

m . 

L 

Given  that  m = x,  it  is  not  a straightforward  task  to  find  the 

ij 

distribution  of  rn  It  will  be  shown  in  Appendix  A that  the  conditional  dis- 
N 

tribution  of  is  Gaussian,  and  that 


E(mN/rnL  = x)  = x-bo-L  + <bN-bL) 


Var  (mN/mL  “ *>  * "l  + "n 


(11-9) 


(11-10) 


if  it  is  assumed  that  only  earthquakes  (as  opposed  to  explosions)  are  consider- 
ed, and  that  the  seismicity  observed  at  LASA  may  be  expressed  by  (II- 8 ).  The 

2 

correction  factor  - b(T  in  ( II - 9 ) arises  from  the  skew  distribution  of  earth- 
quake  magnitudes,  and  is  typically  around  0.  1 units. 

From  the  above  considerations,  it  is  easy  to  show  that  the 
NORSAR  detection  probability  in  terms  of  LASA  m^  values  becomes 
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The  NORSAR  detection  curve  in  terms  of  LASA  magnitudes  can 


easily  be  estimated  directly,  by  selecting  reference  events  reported  by  LASA 
and  verifying  whether  or  not  a NORSAR  detection  occurred.  (Ringdal  and 
Whitelaw,  1973a,  1973b.) 

d.  Discussion 

When  considering  the  three  ca'es  discussed  above,  it  is  evident 
that  we  would  desire  to  estimate  the  NORSAR  detection  capability  in  terms  of 
a "true"  magnitude  as  discussed  under  a. 

A comparison  of  the  detection  curves  obtained  by  the  three 
methods  discussed  here  is  presented  in  Figure  II-2  for  a typical  situation.  It 
is  assumed  that  both  LASA  and  NORSAR  are  unbiased,  that  M T = 3.  70,  (r^= 

0.  15,  and  that  <r^  = = o.  25.  A value  of  b = 2.  0 has  been  used  (this 

corresponds  to  a slope  of  0.9  in  the  base  10  seismicity  curve). 

It  is  seen  that  the  50  percent  detection  threshold  has  a slight 

positive  bias  for  the  curve  based  on  LASA  m,  . The  90  percent  threshold  in 

I)  1 

terms  of  the  NORSAR  curve  is  significantly  lower  than  the  corresponding 

threshold  based  on  LASA  m values  (by  about  0.  4 m units  in  this  case). 

b b 

while  the  "true  ' 90  percent  threshold  is  about  midway  between  these  two  values. 

Although  we  will  give  a more  detailed  discussion  of  how  to  esti- 
mate detection  curves  in  Sections  III  and  IV  of  this  report,  it  seems  appro- 
priate to  give  an  example  of  the  importance  of  the  above  considerations. 
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IRingdai  and  Whitelaw  (1973b)  estimate  the  "optimum"  90  percent  incremental 
detection  threshold  for  the  NORSAR  SP  array  to  be  slightly  below  mb  = 4.  3 
for  the  Kuriles- Kamchatka  region  (essentially  in  terms  of  LASA  magnitudes). 
Yet  Dungum  and  Husebye  (1974)  report  that  the  operational  NORSAR  90  per- 
cent incremental  detection  threshold  for  the  same  region  is  about  = 4.0 
in  terms  of  NORSAR  magnitudes.  NORSAR  and  LASA  can  be  considered 
mutually  unbiased  for  th'S  region. 

This  apparent  contradiction  stems  from  the  fact  that  two  dif- 
ferent detection  curves  are  estimated  in  the  two  casts,  (method  b and  c 
respectively).  By  the  preceding  considerations,  we  can  relate  both  estimates 
to  the  "true"  detection  curve,  since  it  is  reasonable  to  assume  that  the  situa- 
tion illustrated  in  Figure  11-2  applies. 

In  this  way,  we  find  that  Bungum  and  Husebye's  estimate  cor- 
responds to  a "true"  90  percent  detection  threshold  of  = 4.1  - 4.  2,  while 
Ringdal  and  Whitelaw's  estimate  corresponds  to  a "true"  mfa  = 4.0  - 4.  1. 

Thus  the  apparent  contradiction  is  removed,  and  the  proper  conclusion  ir  that 
the  NORSAR  array  has  an  operational  capability  for  the  Kuriles-Kamchatka 
region  which  is  very  close  to  the  optimum  capability  that  can  be  obtained  us- 
ing the  beamforming/ filtering  detection  algorithm. 

3.  Applications  to  Surface  Wave  Detectability 

The  Gaussian  model  described  previously  may  also  be  applied 

to  the  detection  probability  of  long  period  surface  waves.  An  added  considera 

tion  in  this  case  is  how  to  compare  estimates  in  terms  of  to  estimates  in 

terms  of  m . This  has  previously  been  discussed  in  papers  by  Harley  and 
b 

Heiting  (1972)  and  Lacoss  (1971);  we  will  in  the  following  briefly  summarize 
the  relevant  considerations.  Our  assumptions  are: 

• Given  an  event  of  true  bodywave  magnitude  m=x,  then  the  true 

surface  wave  magnitude  M is  given  by 


. i 


. J 
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M = f(x)  + r 


(Ii-13) 


where  £ is  a deterministic  function  and  r is  sampled  from 

2 

a normal  distribution  of  zero  mean  and  variance  tr  . 

The  NORSAR  surface  wave  detection  threshold  magnitude  M 

2 1 

i3  normally  distributed;  ~ 0"^,^.). 

Given  that  the  true  surface  wave  magnitude  M = y , then  the 

2 

NORSAR  surface  wave  magnitude  is  N(v  + b^^ , <r 

M and  M„,  are  independent. 

T N 


a.  NORSAR  LP  Detection  Carve  in  Terms  of  True  Surface 
Wave  Magnitude. 

Given  that  M = y , it  is  found  in  a way  similar  to  Example  2a 


that  the  detection  probability  becomes: 


P(  Detect/ M=y)  = <t> 


y - U 


+ b \ 
\A  T MN 

*3  / 


(H-14) 


where 


2 2,2 

(r  = (T  + (T 

3 MN  MT 


(11-15) 


b.  NORSAR  LP  Detectijn  Curve  in  Terms  of  True  Bodywave 
Magnitude. 

2 

Given  that  m=x  , it  follows  that  M is  N(f(x),  o^).  It  is  then 
found  that  the  conditional  distribution  of  M^  is  Gaussian,  and  that 


E(MN/m=x)  = Kx>  + bMN 


(11-16) 


11-11 


Var  (M^/  m=x) 


(11-17) 


2 2 
= ’MN*  "r 


This  implies  that  the  NORSAR  LP  detection  probability  in  terms  of  true  body- 
wave  magnitude  m becomes: 


where 


P(  Detect/  m = x) 


= <t> 


/fW-MMT*hMN 

\ "4 


2 
( T 


r 


(11-18) 


(11-19) 


By  comparing  expression  (11-18)  to  (11-14)  it  is  seen  that  it  is 
not  correct  to  transfer  between  the  two  detection  curves  si  mply  by  substituting 
y = f(x);  one  must  also  take  the  inherent  scatter  ir.  the  relationship 

into  account. 

In  its  simplest  form  the  relationship  between  and  m^  is 

approximated  by  a linear  function  f , such  as  the  Gutenbe  rg-Richter  (1956) 
relation: 

M s 1,59  mL  - 3.97  (11-20) 

s b 

In  this  linear  case,  the  detection  curve  as  a function  of  becomes  Gaussian. 
However,  if  f is  assumed  to  be  a non-linear  function  (Evernden  et.  al.  , (1971) 
and  Tsai  (1972)),  this  no  longer  holds  true.  In  order  to  obtain  a Gaussian 
curve,  in  this  case,  it  is  necessary  to  plot  detectability  as  a function  of  f(m^). 
This  of  course  requires  that  f is  known  a priori. 

Similar  considerations  as  given  above  apply  to  detection  curves 
in  terms  of  magnitude  estimates  that  are  not  necessarily  "true"  (such  as  LASA 
magnitudes).  We  will  not  elaborate  this  point  any  further  here. 
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D.  A NETWORK  OF  STATIONS 

It  is  not  a straightforward  task  to  apply  the  Gaussian  model 
to  jeismic  detection  probability  using  a world-wide  network.  Clearly  the  de- 
tection probability  is  a function  of  the  network  detection  algorithm;  the  most 
common  algorithm  is  to  declare  an  event  if  at  least  M individual  stations 
(out  of  a total  of  N)  show  cor -esponding  detections. 

Under  the  assumptions  presented  earlier,  let  us  assume  that 
the  probability  of  detection  for  station  number  i , given  true  event  magni- 
tude m " x , is 


P^(Detect/m=x)  = <t> 


R) 


i = 1 , 2 , . . . N 


(11-21) 


If  we  set  the  network  detection  requirement  to  M=l,  and  assume  that  differ- 
ent stations  detect  the  event  independently,  we  obtain  the  following  network 
detection  probability. 


P(Network  Detection/  m=x)  = 1 


(11-22) 


Similar  expressions  can  readily  be  obtained  for  alternative 
values  of  M , using  binomial  expansion. 

Clearly,  the  resulting  detection  curve  is  not  a cumulative 
Gaussian  distribution.  An  interesting  question  is  how  well  it  can  be  approxi- 
mated by  a Gaussian  curve,  and  some  insight  into  this  is  provided  in  Figure 
11-3.  This  figure  shows  the  detection  curve  for  a network  of  N stations 
(N  = 1,  2,  10,  100)  according  to  (11-22),  assuming  that  ail  M.  = 4.  5,  <r,=  0.4 
i = 1,  2,  . . . N.  The  dashed  curves  are  Gaussian  distribution  functions  with 
the  same  50  percent  and  90  percent  detection  levels  as  the  theoretical  curves. 
It  is  seen  that  the  Gaussian  distributions  give  quite  good  approximations  over 
limited  magnitude  ranges,  but  do  not  seem  to  fit  the  entire  curve. 
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Part  of  the  reason  that  the  fit  is  not  good  at  low  magnitudes  is 
the  lack  of  symmetry  around  the  50  percent  level  for  the  theoretical  curves. 
Actually,  the  case  of  M=1  gives  the  least  symmetric  detection  curve;  so 
that  the  Gaussian  approximations  will  be  at  least  as  good  if  multi-station  de- 
tection is  required.  The  best  approximation  may  be  expected  for  M = N/2, 
when  the  theoretical  curve  becomes  symmetric  around  the  50  percent  level. 

In  practical  situations,  the  network  detection  probability  func- 
tion is  further  complicated  by  variable  detection  capability  of  individual  sta- 
tions and  distance  factor  variations.  Furthermore,  occasional  outages  of  one 
or  more  stations  in  the  network  may  cause  loss  of  detections.  It  is  clearly 
difficult  to  give  general  statements  as  to  the  applicability  of  the  Gaussian 
model  to  practical  network  situations.  However,  for  a reasonably  stable, 
homogeneous  network,  it  might  be  anticipated  that  the  Gaussian  model  can  be 
used  as  an  approximation  to  the  detection  curve  over  limited  magnitude  ranges. 
An  example  of  this  will  be  studied  in  some  detail  in  Section  IV. 

Finally,  it  might  be  observed  that  the  network  situation  pre- 
sents additional  problems  if  it  is  desired  to  find  the  network  detection  pro- 
babilities in  terms  of  the  network's  own  magnitudes.  (Similarly  to  Example 
2b  in  this  subsection.  ) This  is  because  network  magnitudes  are  usually  com- 
puted by  averaging  the  magnitudes  of  the  detecting  stations.  Thus  the  network 
magnitudes  will  be  biased  high  for  events  of  low  (true)  magnitude,  that  are  not 
detected  by  all  stations.  Better  methods  of  network  magnitude  determination 
(Herrin  and  Tucker,  1972)  would  be  required  to  overcome  this  problem. 
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SECTION  III 

MAXIMUM  LIKELIHOOD  ESTIMATION 


This  section  presents  a detailed  description  of  the  maximum 
likelihood  method  for  direct  detection  threshold  estimation.  Asymptotic  con- 
fidence limits  are  computed  for  the  estimated  parameters,  and  the  validity  of 
using  these  as  approximations  in  practical  cases  is  investigated  by  performing 
computer  simulation.  A brief  description  of  the  indirect  estimation  technique 
developed  by  Lacoss  and  Kelly  (1969)  is  also  included. 


A.  DIRECT  ESTIMATION  METHOD  - GENERAL  DESCRIPTION 

The  purpose  of  the  estimation  process  described  below  is  to 
estimate  the  "detection  curve"  for  a seismic  station  or  network;  i.  e.  , the 
probability  of  detection  as  a function  of  event  magnitude. 


The  basic  assumption  is  that  the  detection  curve  belongs  to 
some  general  class  of  functions,  and  can  be  completely  characterized  by  the 
values  of  a set  of  parameters.  In  particular,  we  will  deal  with  the  Gaussian 
model  which  was  described  in  Section  II.  We  recall  that,  in  this  model,  the 
detection  curve  is  of  the  general  form; 


P(m) 


(Zmr 


(t  -d) 


2 (T 


dt 


(HI-  1 ) 


Thus,  in  the  Gaussian  case,  the  station  detection  potential  is  cnaracterized 
by  the  actual  values  of  the  parameters  fJ.  and  it  . The  problem  therefore 
is  to  estimate  these  two  parameters. 
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The  general  procedure  in  estimating  the  parameters  of  the  de- 
tection curve  for  a seismic  station  is  as  follows: 

• Obtain  a reference  set  of  randomly  selected  events  of  various 

magnitudes  (as  reported  by  the  reference  source). 


For  each  event  in  the  reference  set,  make  a decision  as  to 
whether  or  not  the  station  has  detected  this  event. 

Establish  the  likelihood  function  for  the  observed  pattern  of 
decisions;  detection  versus  no  detection,  using  the  general 
form  of  the  detection  curve. 

Find  the  set  of  parameter  values  of  the  detection  curve  that 
maximizes  the  likelihood  function. 

To  establish  a formal  model,  we  will  assume  that  the  reference 

data  base  consists  of  n seismic  events  of  magnitudes  m ,...m  , respectively. 

1 n 

These  magnitude  values  are  as  reported  by  the  reference  source,  and  must 
therefore  be  considered  as  statistical  estimates  of  the  true  event  magnitudes 
(see  Section  II).  To  simplify  the  notation,  we  will  assume  (with  no  loss  of 
generality)  that  all  the  m^  values  are  different.  (Note  that  m is  a continu- 
ous variable.)  We  thus  arrive  at  a statistical  test  situation,  where  n inde- 
pendent tests  are  carried  out;  for  each  test  the  probability  of  success  (i.  e.  , 
detection)  is  specified  by  equation  (III- 1 ) , for  the  proper  value  of  m . 

Let  x^  = 1 if  the  station  detects  the  i^  event,  i = 1,2,...  n; 

x.  = 0 otherwise.  The  probability  of  a particular  combination  x.  . . . , x 
i In 

occurring  is 

i 

n . 

A (x  . . . x ) = TT  P(m.)Xi  • (1  - P(m.))  " Xi  (III-2) 

i n . , l l 

1=1 

For  a given  outcome  of  this  experiment,  x = a x = a , the  maximum 
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Ull-3). 

For  information  about  the  general  properties  of  maximum 
likelihood  estimators,  we  refer  to  Cramer  (1945).  Before  discussing  further 
the  mathematical  aspects  of  the  estimation  method,  we  find  it  appropriate  to 
comment  briefly  on  some  of  our  assumptions. 

The  most  important  assertion  in  the  above  estimation  method 
(apart  from  the  validity  of  the  general  form  of  the  detection  curve)  is  the  ran- 
domness criterion  in  the  selection  of  reference  events.  The  essential  point 
here  is  to  select  events  independently  of  the  station  for  which  we  want  to  esti- 
mate the  detection  curve.  In  this  way,  the  events  in  the  reference  set,  for 
any  given  magnitude,  will  represent  a randomly  chosen  subset  of  the  total 
number  of  everts  occurring.  Thus  the  percent  detected  will  be  an  unbiased 
estimate  of  the  percentage  that  the  station  would  detect  of  the  whole  event 
population  for  each  magnitude. 

It  is  clear  that  the  procedure  described  here  will  estimate  the 
detection  curve  of  a seismic  station  in  terms  of  magnitudes  from  an  indepen- 
dent station  or  network.  Thus  the  uncertainties  in  these  magnitude  estimates 
affect  the  resulting  detection  curve  as  described  in  Section  II. 

Further,  it  should  be  noted  that  the  actual  decision;  detection 
versus  no  detection  can  be  made  in  several  ways,  such  as 

• Check  to  see  if  the  reference  events  have  been  reported  in  the 

routine  seismic  bulletin  from  the  station  in  question  (thus  esti- 
mating "Operational  detection  threshold"). 
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• Inspect  visually  the  station  waveforms  for  each  reference 
event  to  verify  detection  or  no  detection  (thus  estimating  the 
"event  verification  threshold"  which  might  be  expected  to  be 
superior  to  the  operational  threshold). 

• Apply  new  detection  algorithms  to  the  station  waveform  for 
each  reference  event,  thus  estimating  the  performance  of  these 
algorithms. 

• For  a seismic  network,  determine  the  detection  threshold  as 
a function  of  the  number  of  individual  station  detections  re- 
quired for  network  detection. 

It  is  important  to  remember  that  the  Gaussian  model  (or  other 
models)  are  to  be  considered  only  as  approximations  to  the  true  station  de- 
tection curve.  Tne  estimation  method  described  here  is  in  a sense  mainly 
a curve-fitting  technique,  and  will  work  best  in  the  magnitude  range  where 
the  largest  number  of  events  are  available.  Inference  about  detection  pro- 
babilities obtained  by  extrapolating  the  curve  to  magnitudes  where  data  are 
scarce  should  therefore  be  avoided.  In  many  cases  it  might  be  advisable  to 
select  only  a specific  magnitude  interval  of  interest,  and  fit  the  Gaussian 
curve  to  the  observed  data  in  this  interval.  Examples  of  this  will  be  given 
later. 

B.  DERIVATION  OF  APPROXIMATE  CONFIDENCE  LIMITS 
1.  Asymototic  Properties  of  the  Estimators 

One  of  the  most  prominent  features  of  maximum  likelihood 

estimators  is  that  they  often  possess  very  desirable  asymptotic  properties. 

/■» 

Under  reasonably  general  conditions,  the  following  may  be  proved  (Cramer, 
1945,  pp  500-504). 
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• The  solution  of  the  likelihood  equation  converges  in  probability 
to  the  true  parameter  value  as  the  rumber  of  observations  in- 
creases. 

• The  maximum  likelihood  estimator  is  asymptotically  efficient. 
(Informally,  an  efficient  estimator  is  one  that  has  lower  var- 
iance than  any  other  unbiased  estimator). 

• The  maximum  likelihood  estimator  is  asymptotically  normally 
distributed. 

In  order  to  verify  that  these  properties  apply  to  our  particular 
case,  we  find  it  convenient  to  regard  the  selection  of  reference  events  as  a 
random  experiment.  Thus  we  assume  that  the  probability  density  function  for 
selecting  a reference  event  of  magnitude  m is  of  some  fixed  form,  s(m)  , 
and  that  individual  selections  are  independent. 

In  this  way  we  can  view  our  direct  estimation  procedure  as 
consisting  of  observing  the  outcomes  (m,x)  of  n independent  experiments, 
each  with  the  likelihood  function: 

\j(m,  x;/i,  (t)  = s(m)*P(m)X»  (l-P(m))1  X (111-4 ) 

The  original  likelihood  function  ( III -2 ) is  then  equivalent  to  a product  of  n 
functions  of  the  form  (111-4);  in  the  sense  that  the  factors  originating  from 
s(m)  do  not  depend  upon  H and  <r  , and  therefore  will  not  influence  the 
maximum  likelihood  estimates. 

As  n— ♦ oo  , we  can  now  apply  the  two-dimensional  form  of  the 
limiting  theorem  in  Cramer  (1945),  in  order  to  show  that  the  asymptotic 
properties  described  above  apply  to  our  case. 
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2.  Computation  of  Asymptotic  Confidence  Limits 

Our  next  step  is  to  find  approximate  values  of  the  variances 

A A 

of  M and  cr  and  their  covariance.  In  view  of  the  preceding  considerations, 
it  is  reasonable  to  compute  the  corresponding  quantities  for  an  unbiased,  ef- 
ficient estimator  of  p and  cr  (usually  known  as  the  Cramer -Rao  bounds), 
and  use  these  values  as  approximations.  Following  Cramer  (1945,  pp  490- 
495)  we  obtain: 


[var  (£)]'*  » (1  - P2(^,a))  . E(  — ) 


(HI-5) 


[var  (o')  ] a (1  - nZ(h,a))  • E ( -d  ^ L ) (IH-6) 


Cov  (fa,  a)  = P(p.a)  • [var(p)  . Var  (a) 


1/2 


< III -7 ) 

In  the  above  formulae,  E denotes  the  statistical  expectation,  with  the  sam- 
ple space  consisting  of  all  possible  outcomes  x,,.  , . x of  detections/no 

1 n 

detections,  each  combination  of  x^'s  having  a probability  of  occurrence  de- 
fined by  equation  (III-2). 

The  correlation  coefficient  P(p,ir)  used  in  the  above  ex- 
pressions can  be  approximated  by 
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P(M,cr)  -E 
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(III-8) 


Because  of  the  simple  form  of  the  likelihood  functions  (III-2) 
and  ( III- 3 ),  it  is  relatively  easy  to  obtain  analytic  expressions  of  the  quantities 


III -6 


f 


defined  in  equations  ( III- 5)  through  (III-8).  By  carrying  out  the  necessary 
computations,  and  in  particular  observing  the  independence  of  the  n test 
situations,  we  find  (setting  P(rn^)  = P^): 


Note  that  expressions  (III- 9 through  III- 11)  are  not  restricted 
to  the  Gaussian  model  (III- 1 ),  and  that  extension  to  classes  of  detection  curves 
with  more  than  two  parameters  obviously  is  possible. 

In  the  Gaussian  situation,  the  function  P is  given  by  (III- 1 ), 
and  we  obtain  after  evaluating  the  partial  derivatives: 
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■ E 

' i=l 

E(J^  . . Z (AA) 
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■ t )‘ • «, 

i=l 


(III- 12) 


Z.  (Ill- 13) 


(III- 14) 
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whr  r n : 


2 2 f 2 

z.  = exp(-(m.-/i)  /a  ) . 27r<r  P(rn.)  • (1  - P(m. 


(Ill-  15) 
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Unce  the  variances  of  and  cr  and  their  covariance  are  known, 
it  becomes  easy  to  find  the  uncertainty  involved  in  estimation  e.  g.  , of  the  90 
percent  incremental  detection  threshold.  In  fact,  according  to  the  Gaussian 
model  (III-I),  ‘he  90  percent  limit  fi  is  given  by 


"90  ' " + 1 • " 


(III- 16) 


where  f = 1.28  is  the  90  percentile  of  the  standard  normal  distribution.  Thus 
we  obtain  an  estimate  anc*  *ts  variance  as  follows; 


A A 

= + f • (T 


(III-  17) 


Var(/i90)  = Var  + 1 * Var  (<r)  + 2f  • Cov  (III- 18 ) 

A 

By  assuming  that  is  approximately  normally  distributed,  we  can  then 

use  ( III  —17)  and  (III-I8)  to  obtain  confidence  limits  on  U 

90 

It  is  in  this  way  possible  to  construct  "confidence  curves"  for 
the  detectability  curve  (III-I)  as  shown  in  an  example  in  Figures  HI- 1 and 
III-2.  Figure  III- 1 gives  a hypothetical  combination  of  events  and  detection 
status.  The  corresponding  maximum  likelihood  detection  curve  and  its  90 
percent  confident  e limits  are  shown  in  Figure  III-2.  In  this  particular  case, 
most  events  in  thi  reference  set  have  magnitudes  greater  than  the  apparent 
value  of  the  50%  detectability  limit  H . Thus,  as  can  be  expected,  the  con- 
fidence in  the  90  percent  detection  estimate  is  greater  than  e.  g.  , the  confi- 
dence in  the  10  percent  estimate. 
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It  is  important  to  interprete  these  confidence  curves  properly. 
First  of  all,  is  should  be  pointed  out  that  the  confidence  linvis  are  derived 
without  reference  to  how  "well"  the  experimental  data  actua’ly  fit  the  Gaussian 
model  (III-  1 ).  The  term  confidence  in  this  case  just  means  U at  if  the  detec- 
tion probabilities  are  according  to  the  model  (111-1),  then  we  can  expect  a 
maximum  likelihood  estimated  based  on,  say,  100  events  to  reflect  the  true 
parameter  value  with  the  given  level  of  confidence.  However,  if,  for  some 
reason,  the  Gaussian  model  (III- 1 ) does  not  appropriately  describe  the  situa- 
tion, our  maximum  likelihood  estimate  as  well  as  its  confidence  limits  will 
of  course  be  meaningless. 

Secondly,  the  confidence  limits  represent  a smoothing  over 
all  data  points.  Thus  there  is  no  reason  to  expect  90  percent  of  all  the  ob- 
served percentages  (marked  as  asterisks)  to  lie  within  the  90  percent  confi- 
dence limits.  This  becomes  more  evident  if  we  decrease  the  size  of  the  mag- 
nitude bins  such  that  only  one  event  corresponds  to  each  magnitude.  In  that 
case  all  "observed"  percentages  will  be  either  0 or  100;  thus  all  correspond- 
ing points  will  lie  outside  the  90  percent  confidence  limits. 

C.  SIMULATION 

The  properties  of  our  maximum  likelihood  estimators  derived 
previously  in  this  section  are  only  asymptotically  valid.  This  means  that 
although  the  given  approximations  can  be  expected  to  work  well  for  large 
reference  event  samples,  nothing  is  said  about  their  validity  when  the  number 
of  events  is  limited.  It  is  clearly  important  to  obtain  some  information  about 
how  well  the  given  expressions  approximate  the  real  distribution  of  the  esti- 
mators in  practical  situations,  and  a convenient  way  to  do  this  is  to  establish 
a simulation  model.  The  simulation  procedure  is  as  follows: 

• Assume  that  the  detection  curve  is  known;  e.  g.  , that  the  Gau- 

ssian model  (111- 1 ) is  valid  and  that  M = M , & = ct  . 

o o 
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Assume  further  that  the  reference  event  set  is  specified;  i.  e. , 


that  the  magnitude  distribution  m,,  . . . m is  known. 

In 


Under  these  assumptions,  we  are  able  to  proceed  as  follows: 


Simulate  the  outcomes  of  100  test  situations,  where  each  sit- 


uation consists  of  the  following  two  steps: 

Creute  a pattern  of  decisions;  detection/ no  detection  for 
the  reference  events,  using  the  predefined  detection  pro- 


babilities in  a random  model. 


A A 

Find  maximum  likelihood  estimates  n and  <r  based 


upon  the  outcome. 


Compute  the  theoretical  approximate  90  percent  confidence 
ellipse  for  the  maximum  likelihood  csf  -nators;  based  upon  the 


equations  (III-5)  through  ( III -8) , with  M = M , rr  = J 

o o 


Compare  the  observed  (simulated)  outcomes  with  the  theoretical 
confidence  ellipse  to  determine  how  well  the  approximations 
work  in  practice. 


The  theoretical  90  percent  confidence  ellipse  corresponding  to 


two  binormally  distributed  estimators  X and  Y with  correlation  P and 


2 

mean  and  variance  (m  , cr  ) and  (m  , a ) respectively,  is  (Cramer,  1945): 

xx  y > 


2 P(x  - m ) (y  - m ) / y - m 

£ L.  + ( 2 


a2(l  -n?)  (III- 1 9) 


where  a is  given  by 


1 - exp  (-  — ) = 0.90 


(III-20) 


III- 1 1 
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In  our  case,  we  can  apply  these  equations  by  substituting  the 
values  given  by  (111-5)  through  (111-8)  and  approximating  the  distribution  of 

A A 

H and  a by  a joint  normal  distribution. 

Figure  111-3  presents  the  results  from  one  such  simulation 
experiment.  The  reference  event  set  contained  n=100  events,  of  a magni- 
tude distribution  identical  to  that  shown  in  Figure  111-1.  Parameter  values 

of  M = 3.  76  and  a = 0.  4 1 we  re  used  as  "true  " parameters  of  the  detection 
° ° 

curve  (111-1).  The  resulting  100  pairs  of  estimates  (At,r r)  are  plotted  in  Fig- 
ure 111-3  together  with  the  90  percent  confidence  ellipse  based  upon  the  theo- 
retical considerations  presented  earlier. 

It  is  seen  that  the  ellipse  reflects  well  the  simulated  distribu- 
tion of  the  estimators  in  this  case.  Thus  we  can  conclude  that  the  "confidence 
curves"  of  Figure  111-2  covering  this  particular  case  provide  reasonably  ac- 
curate indications  of  the  uncertainties  in  the  estimated  parameter  values. 

A second  simulation  case  is  presented  in  Figures  111-4  through 
111-6.  Again,  a reference  event  set  is  given,  this  time  with  n=20  events,  of 
a magnitude  distribution  as  shown  in  Figure  111-4.  A hypothetical  combination 
of  detection/ no  detection  decisions  gave  estimates  of  the  detection  curve  as 
shown  in  Figure  111-5.  By  simulating  100  cases  with  A^q=  4.  10,  ir^=  0.  39, 
we  obtained  the  results  presented  in  Figure  II1-6.  The  90  percent  confidence 
ellipse  is  seen  to  contain  85  of  the  observed  values;  thus  reflecting  reason- 
ably well  the  uncertainties  of  the  estimators.  However,  the  distribution  of 
points  is  clearly  not  symmetric,  and  for  eight  test  cases  (marked  as  arrows) 
a very  large  estimate  of  rr  (greater  than  I.  0)  was  found.  This  indicates  a 
lack  of  stability  in  the  estimation  process  which  we  attribute  to  the  low  num- 
ber of  evems  (n=20)  in  the  reference  data  base. 


In  conclusion,  it  is  clear  that  no  general  statement  can  be  madi 
as  to  the  minimum  number  of  evmts  n that  is  required  for  our  approximations 
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to  be  reasonable.  This  is  because  the  performance  of  the  approximations 
depend  on  the  event  magnitude  distribution  relative  to  the  detection  curve 
as  well  as  on  the  number  n . However,  it  may  be  postulated  that,  given  a 
situation  with  at  least  100  reference  events,  our  approximations  may  be 
used  with  confidence,  provided  that  the  magnitude  range  and  the  distribution 
of  detection/no  detection  decisions  are  reasonable.  For  smaller  event  popu- 
lations, the  approximate  confidence  limits  must  be  used  only  with  great  cau- 
tion, and  with  the  understanding  that  the  estimation  process  may  occasionally 
lead  to  results  with  a very  large  error  margin. 


D.  INDIRECT  ESTIMATION  METHOD 

Indirect  estimation  of  detection  thresholds  is,  in  our  termin- 
ology, a two-step  process: 

• The  first  step  is  to  estimate  the  seismicity  of  a certain  region 
over  a certain  time  period,  by  observing  the  number  of  events 
detected  by  a seismic  system  as  a function  of  magnitude. 

• The  second  step  is  to  use  the  estimated  seismicity  curve  toget- 
her with  the  observed  data  to  estimate  detection  thresholds. 

In  practice,  the  seismicity  is  estimated  by  fitting  the  well 
known  frequency  magnitude  distribution  (111-21)  to  the  data: 


N(m) 


a-bm 

e 


(111-21) 


where  N(m)  is  the  number  of  events  of  magnitude  exceeding  m , and  a 
and  b are  parameters  that  characterize  the  seismicity.  (The  corresponding 
base  10  parameters  are  a'  = a/ log  10  and  b'  = b/log  10  (natural  logarithms) 
respectively.  ) 


111-16 


The  actual  curve  fitting  can  be  done  by  a least  squares  tech- 
nique (as  applied  by  liungum  and  Husebye,  1974)  or  by  a maximum  likelihood 
procedure.  We  will  prefer  the  second  approach  in  the  following;  more  spec- 
ifically, we  will  apply  a method  developed  by  Lacoss  and  Kelly  (1969)  for 
joint  determination  of  seismicity  and  detection  parameters. 

1.  The  Likelihood  Function 


Our  assumptions  are  as  follows: 


The  number  of  detected  earthquakes,  K,  from  a specific  re- 


gion and  their  magnitudes  rn  , . . . , m (listed  in  non- 

1 k 


decreasing  order)  have  been  observed  over  a certain  time  in- 
terval by  a seismic  station  or  network. 


The  seismicity  corresponding  to  the  above  region  and  time 
period  is  given  by  (111-21). 


For  given  values  of  a and  b , the  probability  p(m)dm  that 
an  earthquake  occvtrs  with  magnitude  between  m and  m+dm 
is  obtained  by  differentiating  (111-21): 


p(m)  dm  = b • ea  ,3m  • dm 


(111-22) 


and  the  occurrence  of  earthquakes  in  non-overlapping  mag.  itude 
intervals  are  statistically  independent. 


The  likelihood  function  can  now  be  derived  by  assuming  that  a 
and  b are  given,  and  that  the  detection  curve  P(m)  is  known.  We  choose 
the  following  informal  approach,  which  can  be  extended  to  a formalized  proof 
if  desired: 


Suppose  that  the  magnitude  axis  is  partitioned  into  small  inter- 
vals of  length  dm  . The  probability  of  detection  an  event  in  the  interval  (m, 

v 

m+dm)  is 


F(m)dm  = P(m)  • p(m)dm 


(111-23) 


The  probability  of  detecting  precisely  events  in  the  magni- 
tude ranges  (m.  , m.  + dm.)  i = 1,  2,  . . . K,  respectively,  is  giv'.,  by: 

L(K,  m , . . . , m^/a,  b,  P)dm^,  . . . dm^  = 

nF(m.)dm.  . n d-F(m  )dm  ) J 2 

1=1  ail  other 

intervals 

The  second  part  of  the  right-hand  side  of  (111-24)  is  an  infinite  product,  and 
we  can  clearly  replace  it  by  a similar  product  representing  all  intervals  on 
the  m axis  without  changing  its  limiting  value.  We  thus  obtain: 

Iim  XI  " F(m.)dm.)  = lim  exp  I - \ F(m.)dm.  I 

J J \ an ; j 7 


= exp 


00 

7 


F(m)dm 


(111-25) 


where  we  have  used  the  approximation  e ar  1 + x for  small  values  of  x.  The 
integral  in  (111-25)  is  seen  to  be  the  total  expected  number  of  detected  events 
N,  given  a,  b,  anc  P. 

«o 


N = / 


b • e&  • P(m)  dm  (III-26) 


It  is,  of  course,  assumed  that  this  integral  converges. 

The  expression  (111-24)  for  the  likelihood  function  thus 

becomes: 


111-18 


L(K,  m , . . . m / a,  b,  P)  = e N • b,e&  ^ • P(m  ) (HI-27) 

1 K i=l 

m < m < . . . < m 
1 2 K 

which  is  identical  to  the  expression  found  by  Lacoss  and  Kelly  (1969),  using  the 

0 

Poisson  distribution. 

We  would  like  to  point  out  that  our  derivation  of  the  likelihood 
function  does  not  assume  that  the  observed  number  of  earthquakes  follows  a 
Poisson  distribution.  This  is  very  important,  since  it  has  repeatedly  been 
observed  that  the  occurrence  of  earthquakes  cannot  be  adequately  represent- 
ed as  a Poisson  process.  Figure  111-7  shows  as  an  example  the  actual  dis- 
tribution of  the  daily  number  of  events  reported  by  PDE  during  the  year  1972, 
and  compares  it  to  the  Poisson  distribution  with  the  same  mean.  The  two 
distributions  clearly  do  not  match. 

However,  it  is  interesting  to  notice  that  when  integrating  equa- 
tion ( III— 2 7 ) in  order  to  obtain  the  marginal  distribution  of  the  numt  er  of  re- 
ported events  K , given  a,  b,  and  P,  we  find  as  a result  the  Poisson  distribu- 
tion, with  parameter  N . This  is  not  a contradiction,  it  only  means  that  if 
(hypothetically)  our  experiment  could  be  repeated  a number  of  times  with 
fixed  seismicity  parameters  a and  b ; and  with  a constant  detection  proba- 
bility function  P , then  the  total  number  of  observed  events  K would  follow 
a Poisson  distribution. 

In  practice,  of  course,  the  seismicity  parameters  are  highly 
variable  as  a function  of  time  (in  particular  the  parameter  a);  this  explains 
why  natural  seismicity  does  not  follow  the  Poisson  law. 

It  follows  from  the  preceding  considerations  (by  appropriate 
choice  of  the  function  P ) that  the  number  of  earthquakes  in  any  given  magni- 
tude range  follows  a Poisson  distribution  in  the  hypothetical  case  of  constant 
seismicity. 


*11-  1 9 


Derivation  of  Estimators 


We  will  from  now  on  assume  that  the  detection  curve  is  Gau- 
ssian, i.  e.  , that  P(m)  is  given  by  (111-1).  This  implies  that  the  expected 
number  of  events  N , is: 


N = e 


2 2 

a-bfi  + b a-  /2 


(111-28) 


Following  Lacoss  and  Kelly  (1969),  we  obtain  the  maximum 

A A 

likelihood  estimates  a and  b of  a and  b respectively  as  follows  (ex- 
pressed as  functions  of  n and  rr  ): 


a 1 / 4<r 

1/b  = — (m  -H)  1 + / 1 + j 

2 V (S.  -«  ) 


(111-29) 


u2  2 
b (T 


a = b • H + log  K - — - 


(111-30) 


where  m is  the  aveiage  of  the  observed  magnitudes.  Substituting  these 
values  back  into  the  likelihood  function  (111-27)  then  yields  a function  of  and 
r r that  can  be  maximized  by  a computer  procedure.  The  values  of  y.  and  <r 
that  maximize  this  function  are  the  desired  maximum  likelihood  estimates 
of  the  parameters  of  the  detection  curve. 

3.  A Simulation  Case 

In  order  to  obtain  some  information  about  the  fluctuations  in 
the  parameter  estimates  obtained  by  the  indirect  maximum  likelihood  method, 
we  carried  out  a simulation  experiment  as  follows: 

Suppose  that  an  experiment  has  resulted  in  estimate  values  of 
a,  b,  [1  and  (T . Keeping  these  values  fixed,  we  can  simulate  the  outcome  of 
a number  of  similar  experiments  (e.  g.  , 100),  using  the  Poisson  distribution 
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as  a random  model  to  obtain  the  number  of  detected  events  in  each  magnitude 
bin.  (This  is  acceptable  since  the  Poisson  assumption  is  valid  when  the  sys- 


tem  parameters  are  fixed.  ) Subsequently  we  can  estimate  the  parameter 
values  based  on  the  outcomes  of  the  experiments,  and  see  how  they  fluctuate 
relative  to  the  original,  "true"  values. 

Results  from  a simulation  experiment  are  presented  (for  the 
parameters  jj.  and  tr  ) in  Figure  111-8.  The  original  parameters  were  - 
3.  91,  (T  = 0.  12,  a'  = 6.  00  (base  10),  b'  = 1.  00  (base  10).  This  corresponds 
to  an  expected  number  of  detections  N = 133.  A total  of  100  cases  were  sim- 
ulated. The  results  indicate  the  following: 

• 90  percent  of  the  estimates  of  [i  are  between  3.  85  and  4.  00. 

• 90  percent  of  the  estimates  of  cr  are  between  0.  08  and  0.  17. 

• 90  percent  of  the  estimates  of  b'  are  between  0.  83  and  1. 25 

(not  shown  on  the  figure). 

It  is  also  observed  that  the  distribution  of  points  appears  to  be 
reasonably  symmetric,  and  that  no  significant  bias  can  be  seen.  As  stated 
by  Lacoss  and  Kelly  (1969),  the  estimators  process  the  desirable  asymptotic 
properties  described  previously  for  the  direct  estimation  case;  thus  we  know 
that  the  performance  of  the  estimates  will  improve  as  N becomes  large.  It 
is  thus  concluded  that  the  method  can  be  applied  with  good  confidence  if  sev- 
eral hundred  events  are  available.  However,  we  stress  again  the  importance 
of  verifying  the  sta*istical  assumptions  before  applying  the  method,  since 
outliers  in  the  event  distribution  can  very  easily  distort  the  resulting  estimates. 
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FIGURE  III-8 

INDIRECT  MAXIMUM  LIKELIHOOD  ESTIMATES  OF  DETECTION  PARAMETERS 

(jl.a)  BASED  UPON  100  SIMULATED  CASES 


SECTION  IV 
DATA  ANALYSIS 


This  section  presents  some  examples  of  detection  threshold 
estimation,  with  emphasis  on  the  direct  maximum  likelihood  method.  In 
addition,  a comparison  is  carried  out  between  the  direct  and  indirect  methods 
of  estimation,  based  upon  data  from  two  earthquake  aftershock  sequences  re- 
corded by  the  NORSAR  and  LASA  arrays. 

A.  DIRECT  ESTIMATION  METHOD 

In  this  subsection  we  give  several  examples  of  application  of 
the  direct  estimation  method  to  seismic  stations  or  networks.  These  exam- 
ples are  all  taken  from  evaluation  reports  published  recently  by  Texas  In- 
struments Incorporated.  The  general  estimation  procedure  has  been  as 

follows: 

• Obtain  a reference  event  set  (typically  up  to  500  events)  from 

one  or  more  of  the  following  sources: 

The  PDE  bulletin  (Preliminary  Determination  of  Epi- 
centers from  the  World-Wide  Seismic  Network). 

The  LASA  weekly  event  summary  issued  by  the  Seismic 
Data  Analysis  Center. 

The  NORSAR  weekly  event  summary  compiled  at  Kjeller, 
Norway. 

The  bulletin  from  the  International  Seismic  Month,  ISM, 
(February  20  - March  19,  1972)  compiled  at  Lincoln 
Laboratories. 


] 
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For  each  event,  compute  the  expected  arrival  time  at  the  sta- 
tion to  be  evaluated,  and  display  the  waveforms  for  this  sta- 
tion for  a time  interval  covering  the  signal  arrival. 

Make  a decision  as  to  whether  or  not  the  station  has  a detec- 
table signal  corresponding  to  this  event.  Usually  this  decision 
is  made  by  the  analyst,  based  upon  visual  inspection  of  the 
filtered  and  beamformed  signal  t races. 

For  selected  regions,  compile  the  number  of  events  detected 
and  not  detected  by  magnitude,  and  apply  the  direct  maximum 
likelihood  estimation  method  to  determine  the  station  detection 
curves  for  those  regions. 

When  applying  the  estimation  method,  it  must  first  be  asserted 
that  the  conditions  for  validity  are  fulfilled.  In  particular,  the  following  points 
should  be  observed: 

• Independence  between  reference  source  and  station  t<  be  eval- 
uated. 

• Validity  of  the  Gau  nodel.  In  particular  the  selected  seis- 

mic region  should  be  sufficiently  limited  in  sixe  so  that  B- 
factor  variation  can  be  ignored.  Also  the  limitations  of  the 
Gaussian  model  for  network  detection  should  be  remembered. 

• Consistent  reference  magnitude  estimates.  This  is  essential 
to  obtain  unbiased  estimates.  While  single  station  magnitude 
variance  can  be  mathematically  compensated  for  (Appendix  A), 
a problem  like  the  PDE  magnitude  bias  for  small  events 
(Herrin  and  Tucker,  1972),  is  essentially  untractable,  and 
therefore  represents  a serious  drawback  in  applying  the  method. 
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• False  alarm  problems.  The  possibility  clearly  exists  that  a 

"refe  ence  event"  is  really  a false  alarm.  Likewise,  it  might 
happen  that  the  analyst  verifies  an  event  by  mistaking  a noise 
burst  for  a signal.  Caution  is  required  to  minimize  bias  cau- 
sed by  this  effect. 

1.  NORSAR  SP  Detection  Threshold 

The  detection  capability  of  the  short  period  NORSAR  array  was 
evaluated  by  Ringdal  and  Whitelaw  (1973b).  Figures  IV-1  and  1V-2  show  the 
results  for  the  two  subregions;  the  Kuriles-Kamchatka  arc  and  the  remainder 
of  the  Eurasian  continent,  respectively. 

Of  special  interest  is  Figure  IV-2,  which  shows  very  few  non- 
detections and  hence  gives  an  estimate  of  very  low  confidence  of  the  50  per- 
cent detection  threshold.  In  fact,  it  might  be  questionable  whether  the  Gau- 
ssian model  can  reasonably  represent  the  detection  curve  in  this  case,  since 
D-factors  vary  significantly  over  the  15°-60°  epicentral  distance  range  for 
this  region.  Therefore  statements  based  on  extrapolation  of  the  Gaussian 
curve  below  the  magnitude  range  of  the  observed  data  should  be  avoided  in 
this  case. 

2.  NORSAR  LP  Detection  Threshold 

The  detection  capability  of  the  NORSAR  LP  array  has  been 
evaluated  by  Laun,  Shen  and  Swindell,  (1973).  The  detection  curve  for  their 
total  earthquake  ensemble  as  a function  of  bodywave  magnitude  is  presented 
in  Figure  IV-3.  It  is  seen  that  the  Gaussian  curve  fits  quite  well  in  spite  of 
the  large  epicentral  distance  variation  involved,  (20-70  degrees). 

As  in  the  case  of  NORSAR  SP  estimation,  it  was  necessary  to 
eliminate  events  that  had  been  reported  with  NORSAR  as  the  primary  source. 
This  is  because  SP  detection  and  LP  detec'ion  at  NORSAR  are  not  independent; 
clearly  both  probabilities  fluctuate  with  the  NOkSAR  microseismic  noise  level. 
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It  was  obviously  necessary  to  eliminate  presumed  explosions 

from  the  even4  i nsemble  in  this  case,  since  their  M - m relationship  if 

s b 

functionally  different  from  that  of  earthquakes. 

3.  VLPE  Network  Detection  Threshold 

Our  final  example  concerns  the  estimation  of  the  network  de- 
tection capability  of  the  Very  Long  Period  Experiment  stations  during  the 
winter  of  1972.  Figure  IV-4  is  taken  from  Lambert  et.  al.  , (1973),  and 
shows  the  Eurasian  LP  detection  capability  for  a network  of  6 stations  as  a 
function  of  bodywave  magnitude. 

In  view  of  the  considerations  from  Section  II,  it  is  reasonable 
to  ask  whether  the  expected  lack  of  symmetry  in  the  theoretical  network  de- 
tection curve  has  adversely  affected  the  threshold  estimates.  To  investigate 
this,  we  performed  a series  of  estimation  procedures  for  subsets  of  the  ori- 
ginal reference  pvent  set,  and  the  results  are  shown  in  Table  IV-1.  It  is 
seen  from  this  table  that  by  ignoring  reference  events  of  low  magnitudes,  the 
50  percent  threshold  estimate  tends  to  increase  slightly,  while  the  90  percent 
estimate  decreases.  This  trend  is  fairly  stable  up  to  a cutoff  point  of  about 

m.  = 4.  5;  if  more  events  are  eliminate'.!,  the  reliability  of  the  estimates  de- 
fa 

creases  sharply  due  to  the  low  numbe.  of  non-detections. 

Thus  we  conclude  that  one  must  show  caution  when  applying 
the  Gaussian  model  to  a network  situation.  If  one  is  primarily  interested  in 
the  90  percent  threshold,  a possible  approach  would  be  to  use  the  Gaussian 
model  to  estimate  the  upper  part  of  the  detection  curve,  by  ignoring  low  mag- 
nitude reference  events.  It  would  seem  reasonable  to  choose  a cutoff  point 
somewhere  around  the  50  percent  detection  threshold  n that  case. 

It  should  be  added  that  theoretical  considerations  in  this  parti- 
cular case  arc  further  aggravated  by  the  failure  ol  stations  within  the  VLPE 
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network  to  remain  operational  daring  the  entire  time  frame.  Thus,  for  a 
significant  part  of  the  time,  only  one  of  the  six  V LPE  stations  was  operation- 
al; while  the  most  typical  ituation  was  two  or  three  stations  operational 
(Lambert,  et  al.  , 1973). 

B.  COMPARISON  OF  THE  DIRECT  AND  INDIRECT  ESTIMATION 
METHODS 

As  was  stated  in  Section  II,  the  direct  and  indirect  estimation 
methods  estimate  different  detection  curves,  so  that  different  results  must 
be  expected.  In  the  following  we  will  try  to  verify  the  considerations  from 
Section  II  by  applying  the  two  methods  to  identical  test  situations  and  com- 
pare the  results. 

Two  large  aftershock  sequences  were  selected  for  this  purpose; 
one  from  South  of  Honshu,  Japan  December  3-20,  1972  and  one  from  Kurile 
Islands / Hokkaido  region  June  17-30,  1973.  We  chose  to  estimate  the  opera- 
tional detection  capability  for  the  NORSAR  array  by  the  two  methods  for  each 
aftershock  sequence,  and  compare  the  results.  Thus  we  considered  an  event 
as  being  detected  by  NORSAR  provided  it  had  been  reported  in  the  NORSAR 
seismic  event  summary  compiled  at  Kjeller,  Norway.  It  should  be  noted  that 
swarm  situations  -annot  be  considered  as  typically  representing  natural  seis- 
mic activity.  Nevertheless,  we  have  found  these  two  examples  useful  to  ill- 
ustrate the  theoretical  considerations  presented  earlier. 

1.  Direct  Estimation 

A direct  estimation  of  the  NORSAR  operational  detection  thres- 
h jld  for  these  two  aftershock  sequences  was  carried  out  by  Ringdal  and  White- 
law  (1973b).  They  used  the  SDAC/  LASA  bulletin  as  a reference,  and  thus 
checked  each  event  reported  b / SDAC/ LASA  to  see  if  NORSAR  had  a corres- 
po  iding  detection.  Time  intervals  when  the  NORSAR  array  was  out  of  opera- 
tion were  not  considered.  Thus  the  total  reference  data  base  was  194  events 
for  the  Honshu  swarm  and  364  events  for  the  Kuriles  swarm  (Table  IV-2). 

Out  of  these,  NORSAR  detected  106  and  284  events,  respectively. 
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TABLE  IV -2 

NORSAR  AND  LASA  EVENT  DEI  ECTION  PERFORMANCE  FOR  EARTHQUAKE 
SWARMS  FROM  SOUTH  HONSHU  (DECEMBER  3-20,  1972) 

AND  THE  KURILE  ISLANDS  (JUNE  17-30,  1973) 


Honshu  Swarm 

Kuriles  Swarm 

1 Distance  from  NORSAR  (degrees) 

78 

70 

Distance  from  LASA  (degrees) 

80 

69 

LASA  total  detected  events 

192 

364 

NORSAR  total  detected  events 

133 

452 

Common  events 

106 

284 
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Figures  IV-5  and  IV-6  show  the  distributions  by  magnitude 


(LASA  values)  of  the  events  detected  and  not  detected,  together  with  the 


resulting  maximum  likelihood  detection  curves.  Two  important  differences 
may  be  observed: 


The  NORSAR  array  has  significantly  better  detectability  for 
the  Kuriles  swarm  than  for  the  Honshu  swarm.  This  reflects 
a much  lower  seismic  noise  level  at  NORSAR  during  June  1973 
than  during  December  i972. 


The  Kuriles  detection  curve  hao  a significantly  larger  spread 
than  the  Honshu  curve  ( <r  = 0.47  to  »r  = 0.31). 


This  significant  difference  in  values  of  nr  is  a very  interesting 
point  by  itself,  and  gives  us  an  opportunity  to  verify  one  of  the  basic  assump- 
tions of  the  model  described  in  Section  II.  According  to  equation  (11-12)  of 
that  section,  the  variance  of  the  detection  curve  is  given  by 


2 2 2 2 

(T  = <r  + (T  + (T 

T N L 


(1V-1) 


where  (r^,  denotes  the  standard  deviation  of  the  "threshold  magnitude"  while 


v and  <t  are  the  standard  deviations  of  the  NORSAR  and  LASA  magnitudes, 
IN  Lj 


respectively,  relative  to  a hypothetical  true  magnitude. 


In  order  to  evaluate  the  individual  terms  of  equation  (IV-1),  we 
first  examined  the  variations  in  seismic  noise  level  at  NORSAR  within  the  time 
period  covering  each  event  swarm.  It  was  found  that  the  noise  level  remained 
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essentially  constant  in  both  cases,  so  that  would  be  relatively  insignificant. 


A value  of  <r  = 0.  15  would  be  about  right  in  both  cases. 


The  next  step  was  to  estimate  <r^  and  For  this  purpose, 


we  selected  randomly  50  events  from  each  aftershock  sequence,  and  plotted 
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FIGURE  IV -6 

DIRECT  MAXIMUM  LIKELIHOOD  ESTIMATION  OF  THE  OPERATIONAL 
NORSAR  SP  DETECTION  CURVE  FOR  AN  AFTERSHOCK  SEQUENCE 

FROM  THE  KURILE  ISLANDS 
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NORSAR  against  LASA  as  shown  in  Figures  IV-7  and  IV-8.  We  actu- 
ally selected  the  first  50  events  reported  by  PDE  of  4.  0 in  each  case; 

for  our  purpose  PDE  can  be  considered  to  be  independent  of  LASA  and  NOR- 
SAR. Events  occurring  when  one  array  was  out  of  operation  were  disregard- 
ed. A total  of  four  events  (marked  as  asterisks)  were  not  detected  by  both 
LA5A  and  NORSAR;  the  "missing"  magnitudes  were  then  set  equal  to  the  ap- 
parent "threshold  magnitude".  The  slight  bias  introduced  by  this  procedure 
should  not  be  significant. 

The  most  striking  observation  from  Figures  IV-7  and  IV-8  is 
the  much  larger  spread  between  NORSAR  and  LASA  magnitudes  in  the  latter 
case,  which  corresponds  to  the  Kuriles  swarm.  The  observed  standard  de- 
viations are  0.25  and  0.40,  respectively.  Because  of  the  randomness  in 
event  selection,  it  is  clear  that,  with  the  notation  in  the  figures: 
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so  that  we  obtain  an  immediate  estimate  of  the  remaining  terms  of  equation 
(IV-1). 

Inserting  the  above  numbers  in  equation  (IV-1)  we  obtain  values 
of  o' = 0.29  for  the  Honshu  swarm  and  ( r=  0.43  for  the  Kuriles  swarm,  which 
compare  well  to  the  respective  values  of  0.31  and  0.47  estimated  by  the  maxi- 
mum li'.elihood  method. 

The  important  conclusion  from  the  preceding  considerations 
is  not  so  much  the  numerical  results,  but  rather  the  fact  that  there  is  a 
direct  connection  between  the  spread  in  magnitude  distribution  and  the  spread 
in  the  station  detection  curve  as  estimated  by  the  direct  method.  It  would 
appear  that  the  large  spread  in  the  Kuriles  magnitude  ensemble  reflects  a 
greater  variation  in  source  mechanisms  for  this  event  set  than  among  the 
Honshu  aftershocks. 


i 
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NORSAR 
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NORSAR  m VERSUS  LASA  m FOR  50  EVENTS  RANDOMLY  SELECTED 
b FROM  THE  SOUTH  HONSHU  AFTERSHOCKS 
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FIGURE  IV -8 

NORSAR  m VERSUS  LASA  m FOR  50  EVENTS  RANDOMLY  SELECTED 
b FROM  THE  KURILE  ISLANDS  AFTERSHOCKS 
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Indirect  Estimation 


The  indirect  maximum  likelihood  method  was  also  applied  to 
the  two  aftershock  sequences  described  previously.  The  number  of  events 
reported  in  the  NORSAR  bulletin  in  each  case  is  shown  as  a function  of  mag- 
nitude in  the  histograms  of  Figures  1V-9  and  1V-10.  These  figures  corres- 
pond to  the  Honshu  swarm  and  the  Kuriles  swarm,  respectively.  Since  NOR- 
SAR did  not  report  a magnitude  for  the  main  shock  of  the  Honshu  swarm  (due 
to  system  saturation),  the  LASA  magnitude  for  this  event  was  used  instead. 

The  two  figures  also  list  the  estimated  seismicity  parameters 
and  the  parameters  of  the  detection  curves.  The  incremental  seismicity 
n(m)  which  is  given  by 

n(m)  = b • e&  (1V-3) 

is  shown  as  a dashed  curve.  A curve  showing  the  (incremental)  expected 
number  of  detections  E(m)  is  also  drawn.  This  curve  is  specified  by 

E(m)  = n(m)  • P(m)  (1V-4) 

where  P(m)  is  the  estimated  detection  probability  for  an  event  of  magnitude 
m.  The  "goodness  of  fit"  of  this  curve  to  the  histogram  is  an  indication  as  to 
how  well  the  theoretical  model  with  the  estimated  parameters  actually  fits  the 
experimental  data.  We  observe  that  the  fit  appears  to  be  quite  good  in  both 
cases. 

Confidence  limits  have  not  been  computed  for  the  indirect  es- 
timates. However,  an  idea  about  the  confidence  of  the  estimates  of  the  Hon- 
shu swarm  can  be  obtained  from  the  simulation  experiment  in  Subsection  111-D, 
which  dealt  with  a very  similar  situation.  The  confidence  in  the  parameter 
values  for  the  Kuriles  swarm  should  be  somewhat  better,  since  the  number  of 
events  is  higher  by  a factor  of  three. 
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Comparison  of  Results 


The  estimated  parameter  values  by  the  two  methods  (H,  a and 
90  percent  threshold  M^)  are  summarized  in  Table  IV-3.  The  following 
observations  may  be  noted: 


The  parameter  <r  estimated  by  the  indirect  method  is  signi- 
ficantly lower  in  both  cases  (0.  12  versus  0.  31  and  0.  14  ver- 
sus 0.47,  respectively).  This  is  in  excellent  agreement  with 
our  considerations  in  Section  II,  equations  ( II- 7 ) and  (11-12), 
which  imply  that  the  indirect  method  estimates  the  variation 
(r^,  in  the  threshold  magnitude  only,  and  does  not  relate  to 

the  signal  variations  <r  and  cr  . 

N L 


As  expected  from  the  above  observation,  the  90  percent  detec- 
tion thresholds  estimated  by  the  indirect  method  are  much  low- 
er than  the  direct  estimates  in  both  cases  (0.25  and  0.4  m, 

b 

units,  respectively).  In  fact,  our  considerations  in  Section  11 
imply  that  the  "true"  90  percent  threshold  is  somewhere  be- 
tween the  direct  and  indirect  estimates. 


• The  50  percent  detection  threshold  H is  about  the  same  in  the 

two  cases,  after  compensating  for  the  bias  between  NORSAR 

and  LASA  magnitudes  for  the  Kuriles  swarm.  We  had  actually 

expected  the  direct  estimates  to  be  higher  by  an  amount  of 

2 

b<r  (Figure  II-2);  i.  e.  , about  0.07  m units  for  the  Honshu 
i -*  b 

swarm  and  0.  13  m units  for  the  Kuriles  swarm.  However, 

b 

there  deviations  are  within  the  uncertainties  of  the  mathema- 
tical model. 


In  conclusion,  we  again  stress  that  the  important  point  is  not 
so  much  the  numbers  involved,  but  rather  the  fact  that  the  direct  and  indirect 
methods  estimate  two  different  detection  curves,  and  that  the  resulting  estimates 
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TABLE  IV -3 


COMPARISON  BETWEEN  DETECTION  PARAMETER  ESTIMATES 
OBTAINED  BY  THE  DIRECT  AND  INDIRECT  METHODS  FOR 
THE  NORSAR  SP  OPERATIONAL  PERFORMANCE  FOR 
TWO  AFTERSHOCK  SEQUENCES 
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of  thref  hold  parameters  must  be  interpreted  accordingly.  In  our  opinion, 
the  faik  re  of  the  indirect  method  to  take  the  signal  variation  into  account  is 
a seriov  s drawback  with  this  estimation  technique. 
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SECTION  V 

SUMMARY  AND  CONCLUSIONS 


The  first  topic  of  this  report  was  to  define  the  detection  curve 
of  a seismic  station  or  network  as  the  probability  of  detection  as  a function  of 
event  magnitude.  The  following  observations  were  made: 

• Under  reasonable  assumptions,  the  detection  curve  of  a single 
station  (or  jeismic  array)  for  a limited  region  can  be  approxi- 
mated by  a cumulative  Gaussian  distribution  function.  In  this 
Gaussian  model,  then,  the  parameters  ^ and  cr  of  the  dis- 
tribution completely  define  the  detection  curve. 

• The  Gaussian  model  does  not  theoretically  apply  to  seismic 
networks,  but  may  still  be  useful  as  an  approximation  to  the 
network  detection  curve  within  limited  magnitude  ranges. 

• A very  important  observation  is  that  the  detection  curve  of  a 
seismic  system  varies  with  the  choice  of  reference  magnitudes. 
Thus  a detection  curve  estimated  from  a station's  own  magni- 
tudes tends  to  give  a significantly  lower  90  oercent  detection 
threshold  than  if  a different  station's  magni’udes  are  chosen 

as  reference. 


A maximum  likelihood  method  for  estimating  the  detection 
capability  of  a seismic  station  or  network  in  terms  of  magnitudes  from  an  in- 
dependent reference  system  was  presented.  The  method  is  based  upon  a 
direct  verification  of  detection  or  no  detection  for  a set  of  reference  events. 
Our  presentation  can  be  summarized  as  follows: 
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• The  likelihood  function  for  the  method  was  derived  and  asym- 
ptotic confidence  limits  for  the  estimated  parameters  were 
computed. 

• A simulation  experiment  showed  that  t .e  asymptotic  confidence 
limits  were  good  indications  of  the  stability  of  the  estimates  in 
a test  case  with  100  reference  events  (of  which  75  were  detect- 
ed). A test  case  with  20  reference  events  (10  detections)  indi- 
cated that  the  method  should  be  used  only  with  caution  for  data 
samples  of  this  size. 

• It  was  emphasized  that  the  estimation  procedure  is  only  as 
valid  as  the  model.  The  method  is  sensitive  to  "bad"  data 
points,  such  as  a large  event  that  is  not  detected  or  a very  low 
magnitude  event  *hat  is  detected.  A careful  data  screening  is 
necessa  y to  eliminate  observations  that  either  violate  the  in- 
dependence requirement  or  h e questionable  reference  mag- 
nitudes. Thus,  as  an  example,  the  lack  of  consistency  in  PDE 
m^  estimates  suggests  that  LASA  and  NORSAR  may  in  many 
cases  be  better  suited  as  reference  systems  than  PDE. 

A brief  description  of  the  indirect  maximum  likelihood  estima- 
tion method  developed  by  Lacoss  and  Kelly  (1969)  was  also  included.  A simu- 
lation experiment  showed  that  this  method  gave  reasonable  stable  estimates 
in  a test  case  with  an  expected  number  of  133  events.  Data  screening  in  this 
case  would  be  easier  than  for  the  direct  estimation,  and  the  major  concern 
would  be  to  make  appropriate  limitations  to  the  seismic  region  considered,  so 
that  the  Gaussian  model  is  valid. 

Finally,  examples  of  applications  were  shown,  with  emphasis  on 
the  direct  method.  For  two  earthquake  aftershock  sequences,  a comparison 
was  carried  out  between  the  direct  and  indirect  estimation  method.  The  result 
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was  found  to  be  in  agreement  with  the  theoretical  considerations  regarding  the 
detection  curves. 

In  conclusion,  it  is  felt  that  maximum  likelihood  estimation  is 
a feasible  approach  to  obtaining  estimates  of  the  detection  thresholds  of  seis- 
mic stations  and  networks.  When  choosing  between  the  direct  and  the  indirect 
methods  of  estimation,  we  observe  that  the  latter  method  has  the  following 
two  major  disadvantages: 

0 The  seismicity  estimates  by  the  indirect  method  are  based  upon 

detections  by  the  station  itself,  and  may  not  always  be  reliable. 
For  example,  suppose  we  want  to  estimate  the  NORSAR  opera- 
tional detection  capability  for  a region  with  poor  beam  coverage. 
The  seismicity  estimates  for  this  region  based  on  NORSAR  de- 
tections will  then  clearly  be  biased  low,  thus  causing  the  in- 
direct method  to  estimate  too  high  detection  probabilities. 

• The  indirect  method  fails  to  take  the  signal  variance  into  ac- 

count when  estimating  detection  thresholds.  Therefore  the  90 
perev  » thresholds  found  by  this  method  will  always  be  signifi- 
cantly lower  than  the  "true"  threshold  when  estimating  station 
detection  capability. 


For  the  above  reasons,  we  feel  that  the  direct  method  of  esti- 
mation is  a superior  approach  to  obtaining  reliable  detection  threshold  esti- 
mates. This  method  has  the  added  advantage  of  giving  easy  visual  control 
of  the  results.  However,  the  direct  method  does  require  that  a good  reference 
network  or  station  be  available.  In  the  hypothetical  case  of  a "perfect"  refer- 
ence network,  the  resulting  estimates  from  the  direct  method  would  represent 
the  "true"  detection  probabilities.  In  practical  situations,  the  variance  of  the 
reference  magnitude  estimates  must  be  considered  when  evaluating  the  results. 
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As  in  ail  applications  of  statistical  estimation  theory,  it  is 
necessary  to  do  a careful  data  screening  prior  to  applying  the  mathematical 
tools.  It  is  important  to  remember  that  the  estimators,  being  random  vari- 
ables, sometimes  will  produce  results  that  deviate  significantly  from  the 
true  parameter  values.  Thus,  a careful  interpretation  of  the  results  is  re- 
quired when  applying  the  techniques  described  in  this  report. 
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APPENDIX  A 


STATISTICAL  PROPERTIES  OF  STATION  MAGNITUDE  DISTRIBUTIONS 


We  will  in  the  following  focus  upon  some  statistic.il  properties 
of  station  magnitudes  that  follow  from  the  assumptions  made  previously  in 
this  report.  We  will  in  general  not  refer  to  the  detection  performance  of  any 
seismic  system  during  these  considerations,  thus,  we  assume  that  for  any 
seismic  event  the  station  magnitude  is  defined  regardless  of  whether  or  not 
a signal  is  actually  detected.  To  be  specific,  we  will  refer  to  LASA  and  NOR- 
SAR  station  magnitudes  as  estimates  of  a hypothetical  "true"  magnitude. 
However,  the  results  are  clearly  valid  lor  other  stations  or  networks  that 
satisfy  the  basic  assumptions. 

Our  two  basic  assumptions,  as  stated  in  Section  1,  are: 

• Given  than  an  event  has  true  magnitude  m=x  , then 

The  distribution  of  the  NORSAR  magnitude  m^  is 

N<x  + bN*  <rN) 

The  distribution  of  the  LASA  magnitude  m^  is 

N(x  + bL,  „ZL) 

Th^  two  random  variables  rn  . and  m.  are  independent. 

N L 

• The  number  of  earthquakes  N(x)  exceeding  a given  true  magni- 
tude m=x  may  be  expressed  as- 


N(x)  = e 


a-bx 


( A-  1 ) 


Distribution  of  Tr.'1  Event  Magnitude  Given  Station  Magnitude 


The  first  topic  is  to  ind  what  can  be  said  about  the  true  magni- 
tude of  an  event,  given  c.  g.  , the  v tlue  of  the  LASA  magnitude  m of  that 


event.  I or  this  purpose,  we  find  it  convenient  to  consider  the  true  magnitude 

m as  being  a random  variable,  with  a certain  probability  distribution.  For 

example,  if  we  consider  all  events  of  true  magnitude  greater  than  an  arbitrary 

value  m , then  the  probability  density  function  w(m)  of  m can  be  derived 
o 

from  ( A-l): 


w(m)  = b • e 


-b(m  - m ) 
o 


( A-2) 


The  probability  density  function  (A-l)  is  often  referred  to  as 

the  a priori  distribution  of  m . Our  purpose  is  to  find  the  distribution  of  m, 

eiven  that  m = x.  This  conditional  distribution  is  called  the  a posteriori 
L 

distribution  of  m , given  m^  , and  we  will  denote  it  as  h(m/m^). 

Following  Cramer  (1945)  page  508,  we  obtain: 


h(m/m  ) = 

la* 


w(m)  • 8jJn‘L/m) 


w(m)  . g (ni  /m)dm 

Lj  I-i 


; m > m ( A-  3 ) 

” o 


where  g (m  /m)  according  to  our  first  basic  assumption  is  a Gaussian 
L. 

probability  density  function:  2 


= (2<)'U 


(mL-(m+bL)) 


(A-4) 


When  inserting  (A-2)  and  (A-4)  in  the  expression  (A-3)  we  first 

observe  that  the  term  in  w(m)  involving  m cancels.  Thus  it  is  evident  that 

o 

the  expression  (A-3)  has  a limiting  value  as  n^-t-eo,  since  the  integral  in 
the  denominator  is  convergent.  It  turns  out  that  this  infinite  integral  can  be 
evaluated,  and  we  obtain  after  some  computation: 


I 


h(m/  m ) 

La 


• c 


( A-5) 
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Thus,  we  have  shown  that,  given  that  an  earthquake  has  a LASA  magnitude 
m = y , then  the  true  magnitude  m will  be  normally  distributed: 


E(m/  m = y) 

La 


( A-6) 


Var  (m/  m = y) 

La 


(A-7) 


This  is  a quite  interesting  result,  which  says  in  effect  that  even 
if  LASA  magnitudes  m were  unbiased  for  any  given  event  (i  e.  , b^  = 0),  the 
expected  value  of  the  true  magnitude,  given  the  value  of  , would  be  differ- 

ent from  m . At  first  this  may  appear  to  be  somewhat  surprising,  but  it 
L 

can  intuitively  be  explained  as  resulting  from  the  skew  magnitude  distribution 
of  earthquakes.  As  an  illustration,  consider  the  following  example,  assuming 
b = 0:  Earthquakes  of  true  m = 4.  0 and  m = 5.  0 have  the  same  probability 

of  being  seen  at  LASA  with  m = 4.5.  However,  since  there  arc  many  more 
earthquakes  of  m = 4.  0 than  of  m = 5.  0 , it  follows  that  among  all  the  earth- 
quakes of  m =4.5;  many  more  will  have  a true  m = 4.  0 than  m = 5.  0. 

La 

Therefore,  given  that  m = 4.5,  the  probabilities  favor  m = 4.0  over  m = 5.  0 

L 

as  the  value  of  the  true  magnitude. 

2.  Distribution  of  the  NORSAR  Magnitude,  Given  the  LASA 

Magnitude 

Closely  related  to  the  topic  described  under  1.  is  the  problem 

of  determining  the  probability  distribution  of  NORSAR  magnitude  m , given 

IN 

that  the  LASA  magnitude  for  an  event  is  m = y . 

La 


A-  3 


In  order  to  find  this  distribution,  we  first  recall  that,  given 
m , then  the  true  magnitude  m has  a probability  distribution  h(m/m  ) 

Li  L 

defined  by  (A-5)-  Furthermore,  for  each  given  m , the  probability  density 
function  gjq(mjv/m)  according  to  our  basic  assumptions  is  given  as: 


<mN  - (m  + bN)) 


;N<mN/rn)l  (2”'rN)’1/' 


2<t 


N 


( A-8) 


It  then  follows  that  the  conditional  distribution  f(m^/m  ) of 


m given  m can  be  expressed  as: 

IN  Li 


OD 

f(mN/n.L)  ■ / h(m/mL)  . dm 


( A-9) 


-00 


After  evaluating  this  integral,  we  find  again  a normal  distribution,  with  ex- 
pectation and  variance  as  follows: 


E(mN/nV  = n’L  + (bN-  V * b<,L 


(A*  10) 


,,r  I'V'l1  = °l  + °N 


(A* I 1 ) 


3.  An  Example  of  Application 

The  expressions  (A- 10)  and  (A- 11)  for  the  conditional  expecta- 
tion and  variance  of  NORSAR  magnitude  m , given  LASA  magnitude  m , 

^ Li 

can  easily  be  illustrated  by  considering  a plot  of  NORSAR  magnitudes  versus 
LASA  magnitudes.  Figure  A-  1 presents  such  a plot  for  the  Kurile  Islands 
aftershock  sequence  June  17-20,  197  3.  This  example  was  also  addressed  in 


A-4 


P 


V. 


Section  IV.  All  earthquakes  with  either  m >4.4  or  m > 4.  2 were  in- 

IN  La 


eluded,  provided  both  arrays  were  operational.  For  events  that  were  not 
detected  by  both  arrays,  the  missing  magnitudes  were  set  equal  to  the  ap- 
parent SO  percent  detection  threshold.  (It  is  important  not  to  delete  such 
events,  since  the  resulting  diagram  then  would  not  reflect  the  true  spread 
of  the  magnitude  distribution.)  The  cutoff  limits  4.4  and  4.2  were  chosen 
somewhat  arbitrarily,  witn  the  main  consideration  being  to  reduce  the  num- 


ber of  non-detected  events. 


According  to  equation  (A- 10),  the  regression  curves  of  Figure 
A-l  will  be  parallell,  straight  lines.  We  obtain  from  (A-10)  and  the  sym- 


metric expression  of  E(m  /m  ) that: 

La  IN 


2 2 2 

E(m  -nijm^l  + E(m^-mN/  mN)  = - 1>(  (T  ^ CfL)=  -btf  (A- 12) 


Similarly,  we  obtain  from  ( A- 1 1 ) that 


St.  Dev(m  / m ) = St.  Dev(m  /m  ) - O 

IN  La  La  IN 


(A-13) 


We  recall  that,  when  discussing  this  same  example  in  Section  IV,  we  obtained 
estimates  of  b'=  0.72  (bast  10)  in  Figure  IV-10,  and  <7=  0.40  in  Figure  IV-8, 


using  equation  (1V-2).  Therefore  the  right  hand  side  of  (A-13)  becomes  0.40, 
and  the  ric,ht  hand  side  of  ( A-  1 2 ) become  s -2.3  • 0.72  • 0.16  = -0.  265.  We 
can  now  verify  numerically  that  (A-12)  and  (A-13)  hold  true  in  this  particular 
case  by  evaluating  the  left  hand  side  of  the  two  expressions  di  ectly  from 
Figure  A-l. 


Doing  this,  we  obtain  the  following  results,  when  averaged 


over  all  magnitudes: 


St.  Dev(m  /m  ) = 0.39 

IN  La 


( A-  1 4) 


St.  Dev(m  /m  ) = 0.38 

.Li  IN 


(A- 1 5) 


E(m  -m  / m ) = 0.01 
N L L 


(A- 16) 


E(m  ^-m^/ m^)  = -0.30 


(A- 17) 


It  is  seen  that  (A-14)  and  (A-15)  compare  well  wi«.h  the  expect- 
ed result  of  0.40,  while  the  sume  of  (A-16)  and  (A-17)  is  -0.29,  which  also 
is  close  to  the  expected  number  of  -0.  265. 

Because  of  lack  of  data,  we  cannot  use  this  example  to  verify 
that  the  regression  curves  actually  are  straight  lines,  but  we  may  still  con- 
clude that  the  numerical  result  found  above  gives  good  support  to  the  theo- 
retical model  presented  earlier. 


4.  The  JLASA  Seismicity  Curve  Compared  to  the  True  Curve 

A natural  question  to  ask  at  this  point  is  how  the  distribution  of 
earthquake  magnitudes  affects  the  seismicity  curve  measured  at  a given  seis- 
mic station.  Specifically,  let  us  assume  that  the  "true1'  seismicity  curve  has 
the  form  N(x)  given  by  (A-l),  and  that  the  LASA  curve  is 


NjJy)  = e 


a*  - b*y 


(A-18) 


where  N (y)  is  the  number  of  earthquakes  occurring  with  LASA  magnitude 
exceeding  y. 

We  can  find  N by  evaluating  the  following  integral: 

J-j 


00  I 00 


NJy)  = / / N'(x)  . g^(m^/m=x)  dx  I dm  ^ (A-l  9) 


A-7 


where  g is  given  by  (A-4).  Informally,  the  inner  integral  is  the  sum  over 

all  x of  the  number  of  events  of  true  magnitude  x times  the  probability  that 

the  LAS  A magnitude  is  m , given  x . 

JL* 

Upon  evaluating  the  integral,  we  obtain: 


a + ~ b2  - b • (y  - b^) 

NJy)  = e 


( A-20) 


Thus,  by  comparing  (A-18)  and  (A-20),  we  have  found 


that 


u2  2 

b a 

a*  ~ a + b • b + — — — 

L ^ 


b*  = b 


( A-2 1 ) 
( A-22 ) 


The  conclusion  is  therefore  that  when  estimating  seismicity 
based  on  only  one  station,  the  estimate  of  b will  be  unbiased  relative  to  the 
true  value,  while  this  is  not  usually  true  for  the  estimate  of  a . This  as- 
sumes, of  course,  that  the  estimation  method  itself  does  not  introduce  an 
additional  bias. 


A-8 


